Self-similar scaling in decaying numerical turbulence 
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Decaying turbulence is studied numerically using as initial condition a random flow whose shell- 
integrated energy spectrum increases with wave number k like k q . Alternatively, initial conditions 
are generated from a driven turbulence simulation by simply stopping the driving. It is known 
that the dependence of the decaying energy spectrum on wave number, time, and viscosity can be 
collapsed onto a unique scaling function that depends only on two parameters. This is confirmed 
using three-dimensional simulations and the dependence of the scaling function on its two arguments 
is determined. 
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I. INTRODUCTION 

According to the classical Kolmogorov theory of 1941 
1] , hydrodynamical turbulence is an example of a system 
that is self-similar, i.e., the velocity pattern is supposed 
to look similar when viewed at different degrees of mag- 
nification 0- Of course, real turbulence is not precisely 
self-similar because of intermittency effects that are re- 
sponsible for anomalous scaling, but for the present pur- 
pose such corrections can be regarded as small. 

Two different self-similarity behaviors have been dis- 
cussed in the literature: inertial range self-similarity and 
infrared asymptotic self-similarity that we shall be con- 
cerned with here. The most famous one is probably the 
inertial range self-similarity. Kolmogorov Q showed that 
the velocity difference between two points, Sv, increases 
with scale I such that (Sv) oc £ h , where h = 1/3; see also 
Ref. Q. In other words, when looking at the velocity at 
a magnified scale, x — > ax, where a is the magnifica- 
tion factor, then velocities will only be similar if they are 
rescaled by a factor a , i.e. v — > a h v. 

However, at sufficiently small scales, viscous dissipa- 
tion always destroys the self-similarity. This implies that 
there will be a modification to an otherwise perfect power 
law behavior of the shell integrated energy spectrum, 
E{k). This modification can be described by a universal 
scaling function ip(k, v), which depends on the kinematic 
viscosity v. Thus, one can write 

E(k, v) = k q tp{k, v) (forced turbulence), (1) 
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where q 



(1 + 2h) follows from the normalization 



fE(k,v)dk = \{v 2 ). For sufficiently large Reynolds 
numbers, the energy spectrum has an inertial range with 
h = 1/3, i.e., q = —5/3. This spectrum cuts off at the 
wave number kd — (e/f 3 ) 1 / 4 , where e is the rate of en- 
ergy input. This dependence on v can be used to simplify 
the scaling function to a function that has only one ar- 
gument, i.e., 



^(k,iy) = f(k/k d ) = f(k^ 4 e- 1 ^). 



(2) 



Here, / is a universal function that depends, in addition 
to k, only on the outer scale determined by the geometry 
of the system. 

We now discuss the infrared asymptotic self-similarity, 
i.e., in the following the scaling exponents h and q ap- 
ply no longer to the inertial range, but to the subinertial 
(infrared) range. In the case of decaying turbulence, the 
scaling function also depends on time, i.e., ?/> = ip(k, i, v). 
Furthermore, if> is not a priori universal in the sense that 
its form may depend on the initial spectrum. The spec- 
trum also becomes time dependent, 

E(k,t,v) — k q vb{k,t,v) (decaying turbulence), (3) 

where q depends on the initial condition. If the ini- 
tial condition is restricted to be turbulent so that, prior 
to turning off the forcing, the energy spectrum satisfies 
Eq. (Q), q is expected to be somewhere between 1 and 4; 
see Refs. HQ. 

In a recent paper 0, Ditlevsen, Jensen, and Olesen 
found that the scaling function tp{k, t, v) reduces to a 
two-parametric dependence, 



iP(k,t,v)=g(kt a ,vt b ), 



(4) 



where g = g(x, y) is a new scaling function that has only 
two arguments, and a and b are exponents that depend 
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only on the slope of the infrared part of the initial spec- 
trum. 

Using data from decaying wind tunnel turbulence 
it was possible to show Q that the energy spectra for 
different times can be collapsed onto a single graph by 
plotting k~ q E(k, t, v) versus kt a . The dependence on the 
viscosity and hence on the second argument y of the 
scaling function g(x,y), has been discarded. This may 
be appropriate in the large Reynolds number limit. 

The purpose of the present paper is to determine, us- 
ing numerical simulations, the dependence of g(x, y) on 
both x and y. In a first step we determine the depen- 
dence on x by keeping y constant. This is accomplished 
by letting v vary in such a way that vt b = y = const. 
This can obviously not easily be done in wind tunnel 
turbulence (although v could in principle be changed by 
varying the temperature). In a simulation, changing v 
is of course quite straightforward. The dependence on y 
is determined by integrating over x and considering the 
decay law of kinetic energy. 

II. SCALING IN DECAYING TURBULENCE 

We first recapitulate the derivation presented in 
Ref. Q. The unforced incompressible Navier-Stokes 
equation 

3v \ 

— + v ■ Vv + -Vp = vW 2 v , (5) 
at p 

where p is pressure, is invariant under the transformation 

x ^ ax, v -> a h v, t -> a^'H, v -> a 1+h v. (6) 

In order that Eq. Q can be satisfied, the exponents a 
and b have to take certain values. These exponents can 
be determined by requiring that kt a and vt remain in- 
variant under the scaling transformation QBj, i.e., 

kt a -* (a- 1 k)[a^ 1 - h h} a = kt a , so a = yZT, (?) 
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2n 
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-1.0 


1/2 





1.00 


1.5 


-1.25 


4/9 


1/9 


1.11 


2 


-1.5 


2/5 


1/5 


1.20 


3 


-2.0 


1/3 


1/3 


1.33 


4 


-2.5 


2/7 


3/7 


1.43 



TABLE I: Dependence of secondary scaling parameters on the 
slope q of the initial spectrum. The significance of 2n will be 
explained in Sec. 1111 Dl 



zero point of t corresponds to a time where the energy 
spectrum would have been a pure power law, 

E{k,0,v)~k q . (10) 

Such a spectrum is obviously singular and would have 
infinite energy. We therefore refer to t = as a virtual 
zero point. Near t = the spectrum can therefore not 
be self-similar. On the other hand, if g(x,y) is not nec- 
essarily finite in the limit y — > 0, the above conclusion 
cannot be made. We return to this in Sec. IIII Dl In the 
following we consider the case where time is sufficiently 
far away from zero. 

The validity of Eqs J2J) and Q has already been con- 
firmed using data from wind tunnel experiments |(| where 
the viscosity is low enough so that the second argument 
in g(x,y), y = vt b , can be neglected. One goal of the 
present paper is to demonstrate, using direct simulations, 
that Eq. Q is also valid in the case where the second 
argument, vt b , cannot be neglected. We do this by im- 
plementing in a numerical simulation a time-dependent 
viscosity, v = v{t), such that vt b = const for a given 
value of the initial power-law exponent q. 

As a first test of the scaling relationship we consider the 
decay of fields with initial power law spectra; see Eq. I|1(J|I . 
We will then also test the scaling laws for initial energy 
spectra that are not power laws (Sec. Ifff Efl . 



vt h -> (a 1+h v)[aS l - h h] b = vt b , so b 



1 



(8) 



Translating this into a dependence on q we have, using 
q = -(1 + 2h) and hence h = -(q + l)/2, 



1 



q 



(9) 



Note that 6 = for q — 1, i.e., for initial energy spectra 
that increase linearly with k. In Tabled we have listed 
the scaling parameters for several values of q. 

The limit t — » is problematic. For q > 1, i.e. b > 0, 
both arguments of g(x,y) vanish. Assuming that g(0,0) 
is finite, we can conclude that for t — > the dependence 
of g(kt a , vt b ) on v and k vanishes. This implies that the 



III. COMPARISON WITH SIMULATIONS 

The Navier-Stokes equations for an isothermal and 
weakly compressible fluid are solved in a box with peri- 
odic boundary conditions. We always adopt initial veloc- 
ity fields such that their Mach number is around 1%, so 
compressibility effects can be neglected. We employ the 
Pencil Code [8| which is a higher-order finite-difference 
code using the 2 iV-RK3 scheme of Williamson Q for time 
stepping. The low Re runs presented in this paper were 
done on relatively coarse grids (64 3 ) while the high Re 
runs had a resolution of 256 3 . For further details and 
recent turbulence simulations using the Pencil Code see 
Ref. We begin by studying the evolution of initial 

velocity fields with power-law spectra. 
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A. Initial power-law spectra 

The initial power-law spectra with arbitrary values of q 
are constructed by first generating in real space a random 
velocity field that is 8 correlated in space. Such a velocity 
field corresponds to a k 2 energy spectrum. In Fourier 
space, the velocity field v(k) is then multiplied by a factor 

In the upper panel of Fig. ^ we show the results of 
a numerical experiment where a power spectrum with 
q = 1 decays under the action of constant viscosity. In 
this case we have b — 0, so y — const. This means that 
the spectra collapse onto a single graph when E(k, t, v) 
is divided by k q (= k in this case, because q = 1) and k 
is multiplied by t a (= t 1 / 2 in this case). This is indeed 
the case, see the lower panel of Fig. ^ 
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FIG. 2: Same as Fig.0 but for q = 2 and constant viscosity, 
v = 10~ 3 . The times shown in the upper panel are t = 0, 3, 
9, 26, and 77. Note that the curves for different times do not 
collapse onto a single graph. 
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z^-of for t < t Ici 

v-tcf (t/ttef) b otherwise, 



(11) 



where i/ Te f = v(t rc {) is a constant reference viscosity. At 
early times, t < i rc f, the initial fields were allowed to 
decay under the action of a constant viscosity v rc { so 
as to avoid having to use an excessively large (or even 
infinite) viscosity. The result is shown in Fig. [31 and the 
spectra for the different times collapse reasonably well 
onto a single graph. 



FIG. 1: Decay of the initial energy spectrum E(k, 0, v) ~ k q 
for q = 1 (upper panel) and the corresponding scaling function 
g q — E(k,t,v) /k q versus x — kt a (lower panel). Note the 
collapse of the rescaled spectra for different times (t — 0, 1, 
2.5, 6, 12.5, and 29). For q — 1 the parameter a = 1/2. 

For all other values of q, the second argument y = vt h 
will not be constant and must depend on t. We therefore 
expect that the spectrum will not collapse onto a single 
graph. This is shown in Fig. [21 where we show the spec- 
tra at different times (upper panel) and the attempt to 
collapse them onto a single graph (lower panel). 

Collapse of the spectra obtained at different times is in 
general not possible unless one makes v time dependent 
in such a way as to keep y — v(t)t h constant in time. In 
the following simulations the viscosity is therefore given 



B. Dependence of g(x,y) on the first argument 

It turns out that for a fixed value of y and different 
values of q the scaling function g(x, y) does not quite col- 
lapse onto a single graph and that, therefore, the curves 
for different values of q are distinct. We indicate this 
by a subscript q and write g q (x,y). However, empir- 
ically it turned out that to a good approximation the 
(/-dependence can be removed by rescaling x by a q de- 
pendent factor, i.e., 

x — > x = x(^q T 4)/5. (12) 

Note that for q = 1 we have x — x. In Fig. 0] we show 
g q (x,y) versus x for fixed value of y and three different 
values of q (=1, 2, and 3). Note that the collapse of the 
three curves is reasonably good. 
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FIG. 3: Same as Fig. H but with v = u(t) given by Eq. JTTJ 
with u rc f = 3x 10 -3 and t TC f — 0.1. Note that now the data 
points collapse reasonably well onto a single graph. In this 
figure the times are the same as in Fig. [5] 
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pendence of v(t). Thus, v{t) is proportional to t bv where 
b v = — (1 — <?!/)/ (3 + q v ), which is analogous to Eq. 10. 
The result is shown in Fig.[S]and we see that the best col- 
lapse is indeed achieved when q = q v . We have checked 
that this agreement holds also for different values of q. 
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FIG. 5: Scatter plots for q — 2 and different values of q v . 
Note that the best collapse is achieved for q v — q. 



D. Dependence of g q (x,y) on the second argument 



FIG. 4: Three sets of g(x, y) curves from decay experiments 
with different values of q but the same value of y (= 2xl0~ 2 ). 
The abscissa has been rescaled according to Eq. H'2t to make 
the curves for q = 1 (triangles), q = 2 (squares), and q = 3 
(crosses) collapse onto a single graph. 



C. Modified time dependence of viscosity 

In order to verify the anticipated scaling behavior fur- 
ther, we determine effective values of q and check whether 
these values are consistent with each other. We begin by 
defining an effective value q v that determines the time de- 



Next we consider the temporal decay law of the kinetic 
energy. This allows us to constrain the dependence of 
g{x, y) on y. As usual, the kinetic energy (per unit mass 
and unit volume) can be found by integration, 

/>oo 

Elan(t,v)= E{k,t,v)dk (t>0). (13) 

Jo 

For a given value of y, where y may still be a function of 
t, Eq. I|13fl can be rewritten as an integral over the first 
argument of the scaling function, x — kt a . This gives 

/>oo 

E kin (t,y)=t- a ^ x"g q (x,y)dx, (14) 
Jo 
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where -E^in still depends on y = y(t). Here we have ig- 
nored the fact that in order for g q (x,y) to be indepen- 
dent of q we should rescale the x coordinate by a factor 
(q + 4)/5; see Eq. Q12[l. However, this only corresponds 
to an overall rescaling of the kinetic energy by a factor 
[(q + 4)/5]~^ l?+1 ' ) and is therefore unimportant. 

It is convenient to isolate the main t dependence, 
-Ekin ~ t~ 2n , where 2n = ail + q), and 



1 



We can therefore write 



Eun(t) ~ t~ In g q [vt»), 



(15) 



(16) 



where g = g(y) is a function that only depends on 
one argument and is obtained by integrating out the in- 
dependence of g(x,y), i.e. 



9g(y) 



x q g q {x,y) dx. 



(17) 



The resulting values of 2n are given in Table[Qfor different 
values of q. Note that the basic t~ 2n decay law has also 
been obtained in Ref. 0, but there it was assumed that 
v is negligibly small and that for large values of k one 
has a Kolmogorov spectrum. The basic relation between 
2n and h or q can also be obtained by assuming that the 
rate of dissipation is proportional to u| / £, where £ is the 
integral scale 0). 

In Fig. El we check that the basic decay law is indeed 
mostly governed by the value of q and not by the value 
of q u . For q v = q = 2 the decay has the expected slope 
2n = 1.2 (middle panel). For q v = 1, the viscosity is 
constant and the decay is accelerated, while for q v = 4, 
the viscosity decreases faster than is necessary for keeping 
y constant, and the decay of -Ekin is now slower than what 
is expected based on the value of q. These results confirm 
that the best agreement is achieved for q v = q. We have 
checked that the same is true for q v = q = 3, for example. 

We now turn to the y dependence of g q (y), which can 
be determined by plotting t 2n Eki n versus y; see Fig. 
Within plotting accuracy the results seem to be inde- 
pendent of the value of q. We can therefore drop in the 
following the subscript q on g q (y). 

The results confirm that for small values of v the time 
dependence of the decay law of kinetic energy is weaker: 
g ~ y -0 - 85 for y < 0.003 compared to g ~ y -18 for larger 
values of y. However, there is as yet no evidence that g 
becomes completely independent of y when y — > 0. 

The above results imply that the energy decay law is 
attenuated by a small correction factor for q ^ 1. Con- 
sider, as an example, the case q = 2. The basic decay 
law is Eki n ~ t^ 1 ' 2 , see Eq. I)16|). For q = 2 we have 
b = 1/5, so for small values of v (assuming y < 0.003) 
the exponent has to be corrected by —0.85/5 = —0.17, 
so that the correct decay law is -Ekin ~ t~ 137 . This is 
indeed confirmed by direct inspection of the data. 
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FIG. 6: Energy decay law for q — 2 and different values of 
q v (solid lines). The different slopes are 2n = 1.5, 1.2, and 
0.9. The dashed lines indicate the slope expected if the decay 
law was governed by the value of q v (slopes 1, 1.2, and 1.43; 
sec Table ID. A gain, the best collapse is achieved for q v — q, 
corresponding to 2n — 1.2. 



The fact that g does not seem to go to a finite value 
in the limit y — > is surprising, because it implies that 
there is a viscosity correction to the basic t~ 2n decay law 
even in the limit of vanishing viscosity. Although the 
data do not necessarily allow such an extrapolation, we 
are not aware of any evidence against a finite viscosity 
correction in the large Reynolds number limit. 



E. Turbulent initial conditions 

The results shown in the previous sections demonstrate 
that the scaling law ijSJl successfully describes the decay 
of kinetic energy in the special case of a flow field with an 
initial k q spectrum. A more realistic initial condition is 
a turbulent velocity field which has an energy spectrum 
that is decreasing with increasing k. Nevertheless, there 
is always a subinertial range where values of q between 1 
and 4 are not uncommon. 

In the following we consider the decay of flow fields 
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FIG. 7: Representation of g q (y) obtained by patching to- 
gether the decay laws from different runs with different values 
of q. At early times the decay is not yet self-similar, so the 
data for these times have been ignored in the plot (the re- 
mains of the initial transients can still be seen as little hooks 
in the beginning of each piece) . Each piece of the decay curve 
has been shifted along the ordinate to connect the different 
pieces with each other. The normalization of g is therefore 
arbitrary. 



that are initially statistically stationary. These initial 
fields are produced by applying a random force within 
a band of wave numbers around kf until the work done 
by the forcing is balanced by dissipation. Relatively high 
resolution (256 3 ) and large values of kf are needed in 
order to get a well-defined subinertial range. 
As explained in Sec. |nj the scaling law 

Eq. @ 

is a direct consequence of the scaling properties of the 
unforced Navier-Stokes equations, Eq. JSJ, and will not 
be valid for a flow driven by a general forcing function. 
The statistically stationary state considered here will not 
necessarily be compatible with the Navier-Stokes equa- 
tions. In the following we assume that t = is the time 
when the forcing is stopped, but we must expect there to 
be some readjustment phase before self-similar scaling is 
possible. 

When viscosity can be considered negligible or when 
it is made time dependent according to Eq. (JTTJ, the 
parameter q can, in principle, be determined by two in- 
dependent methods. It can be found by determining the 
spectral slope in the subinertial range or by fitting the en- 
ergy decay to Eq. , as done in Ref. [7j . In addition, of 
course, q could be found completely empirically by trying 
different values until the collapse is best. Unfortunately 
the first approach is difficult since one has to have large 
scale separation between the size of the box and the in- 
tegral (or forcing) scale in order to be able to resolve the 
infrared region. This requires very large resolution. In 
addition, the infrared limit obtained from simulations is 
not very accurate for small k, because only a few modes 
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FIG. 8: Upper panel: initial energy spectrum (solid line) 
together with subsequent energy spectra (dotted lines) ob- 
tained from driven turbulence simulation forced at wave num- 
ber kf = 10. The spectral slope for small wave numbers 
corresponds to q = 1.5. Middle panel: attempt to collapse 
the spectra on a single graph which fails at early times. 
Lower panel: decay of kinetic energy corresponding to a slope 
2n = 10/9 ~ 1.11, confirming the q — 1.5 scaling. 



contribute to the shell-integrated spectrum. 

In Fig.|Hlwe show the result for a turbulence simulation 
that was driven at ki = 10 and the forcing was turned 
off at t = 0. Note that the collapse is relatively poor at 
early times. The collapse improves significantly when the 
turbulence is driven at kf = 30; see Fig. 

The reason for the collapse being much better in the 
case of larger kf is probably related to the facts that the 
local turnover time ~ (u Ima £) is shorter. Thus, self- 
similarity can probably commence much earlier. 

Finally, we note that in our simulations the subinertial 
range slope is q = 1.5 both for large and small values 
of kf. We are not aware of a theoretical explanation for 
this slope, but it is probably related to finite size effects. 
By contrast, in an infinite domain the slope is expected 
to be q — 4 (or = 2), which could be motivated if the 
Loitsyansky (or Saffman) integral were independent of 
time. In that case one would have 2n — 10/7 (or 2n — 
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FIG. 9: Same as Fig. 03 but for k { = 30. Note that the 
collapse is much better - even at earlier times - than for the 
case with k f — 10 (middle panel). 

6/5). 

F. Conclusion 

The results presented above have shown that decay- 
ing hydrodynamic turbulence can be characterized by 
a two-parametric scaling function and that this func- 
tion may well be universal and independent of the ini- 
tial slope q of the spectrum in the infrared limit, i.e. 
in the subinertial range. Although the basic scaling 
behavior has already been confirmed earlier using 
data from wind tunnel turbulence |6j, it was not pos- 
sible to determine the infrared scaling properties of the 
energy spectrum. Indeed, for the smallest wave num- 
bers available from the wind tunnel data the spectrum 
was still a decreasing function of wave number k. This 
is because wind tunnel measurements only allow access 
to one-dimensional spectra which are always monotoni- 



cally decaying. This property follows from the fact that 
for isotropic turbulence the three-dimensional spectrum 
E(k) is related to the one-dimensional spectrum Eiv(k) 
via E(k) — —kdEm/dk, and since E(k) > it fol- 
lows that the one-dimensional spectrum can never in- 
crease with k; see Eq. (7) of Ref. . This is also true 
for the longitudinal and transversal power spectra sepa- 
rately; see, for example, Fig. 6.11 of Ref. Whether or 
not the proper subinertial range of the three-dimensional 
spectrum can be determined from wind tunnel experi- 
ments is unclear. It is therefore important that simula- 
tions can now demonstrate explicitly that the slope of 
the subinertial range spectrum is linked to the scaling 
law derived in Ref. 0]. 

There are obvious extensions of this work to the case 
of decaying magnetohydrodynamics turbulence. Simi- 
lar scaling properties also apply to the magnetic case 
EH El- but the detailed functional dependence of the 
corresponding two-parametric scaling function has not 
yet been fully determined, although partial results do al- 
ready exist. In particular, for the case of helical initial 
fields the combined dependence on wavenumber and time 
has been studied in Ref. [ljj , and resistive corrections to 
the decay law have been investigated in Ref. [lq. This 
work generalizes earlier findings that in the helical case 
the magnetic energy can decay as slowly as ~ i -1 / 2 [l6j . 
while in the nonhelical case the decay is g enerally faster 
and similar to the hydrodynamic case \17\ - 

A general difficulty with the self-similarity approach 
is the uncertainty regarding the zero point of t. There 
is apparently no unique way of determining this time 
a priori. An a priori choice of the zero point of t is 
however necessary if v is allowed to be a function of time. 
Although the uncertainty regarding the zero point of t 
becomes less influential at later times, it is not normally 
possible to revise the zero point of t afterwards, unless 
one is prepared to run an entirely new simulation. 

Once the initial startup phase is over and the decay 
has become self similar, one is however able to determine 
in full detail the exact form of the two-parametric scaling 
law. Our current work can only be preliminary, because 
it remains to be checked how general and perhaps even 
universal the g(x,y) function really is. If its generality 
is established, it could become a powerful analysis tool 
for making predictions about the decay of kinetic energy. 
This applies in particular to the function g(y), which 
plays the role of a viscosity dependent correction function 
for the decay law of the kinetic energy. 
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